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Abstract 

We propose a time discretization scheme for a class of ordinary differential equa- 
tions arising in simulations of fluid/particle flows. The scheme is intended to work 
robustly in the lubrication regime when the distance between two particles immersed 
in the fluid or between a particle and the wall tends to zero. The idea consists 
in introducing a small threshold for the particle-wall distance below which the real 
trajectory of the particle is replaced by an approximated one where the distance is 
kept equal to the threshold value. The error of this approximation is estimated both 
theoretically and by numerical experiments. Our time marching scheme can be easily 
incorporated into a full simulation method where the velocity of the fluid is obtained 
by a numerical solution to Stokes or Navier-Stokes equations. We also provide a 
derivation of the asymptotic expansion for the lubrication force (used in our numer- 
ical experiments) acting on a disk immersed in a Newtonian fluid and approaching 
the wall. The method of this derivation is new and can be easily adapted to other 
cases. 



1 Introduction 

One of the challenges for fluid/particle flow simulations is to provide an accurate resolution 
of the lubrication regime when the distance between two particles immersed in the fluid 
or between a particle and the wall becomes very small. Taking aside the problems related 
to the discretization in space (extremly high gradients of the velocity in the narrow gap 
between the particle and the wall, for example), we focus our attention in this article on 
the discretization in time. The lubrication regime is characterized by very high magnitude 
of the drag force (which can be referred to as the lubrication force in this case). Very 
small step sizes in time should be thus employed in order to obtain a physically acceptable 
solution. Our idea is to prohibit the particle from approaching too closely the wall during a 
simulation. We shall thus choose a threshold q s for the distance q between the particle and 
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the wall and replace the true trajectory of the particle by an approximated one, in which 
the distance q is kept equal to q s until an eventual rebound of the particle from the wall. 
The moment of rebound is predicted using an auxiliary quantity (a crude approximation 
of the velocity) that is computed all along the period of time when the particle is stuck 
at the distance q s . This approach reminds the gluey particle model of [HI [12] where q s 
is set to zero and the limit of vanishing viscosity is considered. However, our motivations 
are quite different from that behind the gluey particle model. This model is intended as 
a simple alternative to the standard governing equations of Navier-Stokes type, eventually 
corrected by taking into account the roughness of the particle surface. On the other hand, 
our approach is to take the standard fluid equations for granted (assuming the particle 
surface to be smooth) and to provide a tool for a robust time discretization of them. It 
means, in particular, that we would need an accurate enough method valid for any given 
value of the viscosity, not necessarily small. Note also, that our threshold q s will typically 
depend on the time step size, so that it is indeed a numerical device and it has no physical 
meaning. 

The plan of the article is as follows: we start by reminding the governing equations in 
a general setting and by explaining in more detail the difficulties related to the simulations 
in the lubrication regime in the next section. Section 3 is the core of the paper. The idea of 
the threshold is rigorously introduced and studied there for the model ordinary differential 
equation representing the essence of the general setting in the simplest case of a circular 
(or spherical) particle approaching the wall. The discussion is held on the continuous level 
in Section 3. The discretization in time is introduced in Section 4 where several implemen- 
tations of our idea are proposed on the discrete level followed by numerical experiments. 
We also include an appendix detailing a derivation of the asymptotic expansion for the 
lubrication force acting on a disk approaching the wall. The method of this derivation is 
new and can be easily adapted to other cases. 

2 Motivations: governing equations for the fluid/particle 
flows and some difficulties arising in their simula- 
tions 

The general setting of this work is a study of the motion of a rigid particle immersed in 
a viscous, incompressible fluid with a particular emphasis on situations when the particle 
approaches the plane. To set the notations, we assume in general that the fluid (with the 
particle inside) fills a fixed domain fl e M. d with d = 2 or 3, while the region occupied by 
the particle Bt C Q varies with time t. We denote the time-dependent fluid domain JF t so 
that = Bt U Tt at any time t. Supposing that the inertial effects are negligible in the 
fluid and the no-slip conditions are valid on the boundaries of Q and B t , the fluid motion 
is governed by the Stokes equations 
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Figure 1: Notations 



' -z/Au + Vp = p/g, in JF t 

V ■ u = 0, in Tt 

u = 0, on dQ 

k u = V + wxr on <9£>, 



(1) 



where u and p are the velocity and the pressure in the fluid, v and pf are the viscosity and 
the density of the fluid, g is the external force, V = V(t) and u — u(t) are the translational 
and angular velocities of the rigid body B t , r = x — G is the vector pointing from the center 
of mass of the particle G to a point x on its boundary. The more realistic Navier-Stokes 
equations may be accomodated into the framework (Tj[|) by including the convective term 
u • Vu into g. 

The fluid exerts a net force F and a torque T on the particle given by 



where D(u) stands for the symmetric gradient of u and n is the unit normal vector on dB t 
directed towards the fluid domain. Note that F and T are indeed functions of only the 
placement of the particle Bt and its translational and angular velocities since the velocity 
u and the pressure p are uniquely determined in the fluid by the parameters Bt, V and u 
as the solution to the Stokes equations ((H). Moreover, the dependence of F and T on V 
and u> is linear. Using these notations we write out the equations of motion of the particle 
as follows 




(2) 




+ uxl t u = T(B t ,V,u), 



(3) 
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where m is the mass of the particle and It is its inertia tensor, expressed in the fixed 
Cartesian frame and thus dependent on time. Equations ([3]) are coupled with the equations 
describing the propagation of the particle, i.e. G = V for the center of mass and r » = wxr,;, 
i — 1, . . . , d, for the vectors r, fixed in the particle. The force F is a sum of the Archimedes 
force due to gravity and of the drag force which is purely hydrodynamic, i.e. obtained 
from ([1]) by setting g = 0. The particularity of the drag is that it tends very rapidly to oo 
when the particle approaches the wall, thus preventing collisions between them. Indeed, 
it has been proved in [8] (2D case), [9] (3D case), that a smooth rigid body embedded in 
a viscous fluid cannot touch the wall in finite time. In the regime of very small distance 
between the particle and the wall, the drag force is also known as the lubrication force and 
it is notoriously difficult to take into account in a numerical simulation. 

It is noteworthy that these considerations are valid only for smooth surfaces; however, 
modelling surface roughness of the wall or particle is a much more delicate issue. Several 
experimental ([IB], [IT] . [TU] . [T5] ) and theoretical ([E], [H]) works propose in this case to 
modify the expression for the lubrication drag force by introducing a shift to the distance 
to the wall, of magnitude strictly lower than the roughness size. This physically means that 
roughness decreases the dissipation in the system, and that the interaction is similar to 
that between equivalent smooth surfaces, located at some intermediate position, between 
the peaks and the valleys of asperities. 

To simulate the motion of the particle governed by equation described above, one should 
be able to solve numerically the Stokes system (CQ) for any given particle position Bt and 
for any given velocities V and u>. This is a formidable task in itself especially because 
the position of the particle is not known a priori but changes with time. We do not make 
precise the choice of the numerical method for the Stokes system. We just assume for the 
moment that the force F and the torque T can be computed for any B t , V, u but this 
computation is in general very expensive. Our primary goal in this article is to devise an 
efficient discretization in time of fl3]) using as few as possible solutions of the Stokes system. 
The simplest idea is to use the following scheme 

~\ t ~\ t k ■ ' 1 

m ^ = F(3fc_ lf V k ,c/) + mg(t fe ), (4) 

1 x \_y = T(iV x , VV fc ), (5) 

= v *. <"> 

yi k y, k 1 

1 . = Y i k XL0\ i = l,...,d. (7) 

We have introduced here the uniform grid in time tk = kAt and have denoted the quantities 
computed at the time step tk by the superscript k. The idea behind the scheme (jl])-((7j) is 
to compute first the velocities by (jl])-© and then to propagate the particle using the last 
available values of the velocities. The equations (jH)-© are thus coupled together and also 
coupled with the solution of the Stokes system on a fixed geometry given by position of the 
particle B th ^ on the previous time step tk-\- The cost of such a computation is normally 
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essentially the same as that of the Stokes system ((T]) with prescribed V and u>. We need thus 
one solution of the Stokes system per time step. This approach was successfully used, for 
example, in [13] in conjunction with a fictitious domain discretization of the Stokes system 
as in However, independently from the discretization in space, one would encounter 
problems when the particle approaches the wall. Indeed, in this lubrication regime the force 
F explodes and thus the scheme (jlJJ-flSj) is no longer valid unless an extremely low value for 
the time step is used that makes the simulation prohibitively expensive. A commonly used 
cure for this problem is to introduce short-range repulsion forces between the particle and 
the wall, as in [6], for example. However, the influence of these (not necessarily realistic) 
forces on the accuracy of a simulation is not well understood. Another simple idea is just 
to stop the particle when it tries to penetrate the wall during a numerical simulation. 
However, it is then not necessarily clear what criterion should be chosen to decide if the 
particle should eventually bounce off the wall and when should it happen. These questions 
have a partial answer in the articles [HJ [T2] on the gluey particle model. It is shown there 
that the particle trajectory satisfies an integro-differential equation in the limit of vanishing 
viscosity, which is easy to discretize in time using moderate time steps and which predicts 
the moment of an eventual rebound from the wall. We pursue a similar idea in this article 
but our aim is to construct an approximated trajectory of the particle in the lubrication 
regime that would be accurate enough for any given value of the viscosity, not necessarily 
small. 

3 A model ordinary differential equation with lubri- 
cation forces 

3.1 The model 

Let us consider the simplest setting of the problem described in the previous section: 
assume Q C M 2 is the half-plane {(x,y), y > 0} and the particle is a disk of radius R. 
Let moreover g(£) = g(t)e2 and assume the particle is at rest at the initial time. The 
^-component of the particle velocity and its angular velocity will then vanish at all time. 
The position of the particle is fully determined by its distance q from the bottom, as in 
Figure [TJ The net force F is the sum of the drag, which is a function of q and V, linear in 
V, and of the Archimedes force: 



where n(q) is the drag coefficient computed by the Stokes equations (JTJ). Denoting the 
y-component of the velocity by v, we are thus led to the following differential equations 



F 



n(q)V + m a g(t)e 2 , with m a = m- p f \B t 
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We are especially interested in the lubrication regime of small q. The asymptotic of n(q) 

3 

when q — > is given by 3y/2nu (f ) 2 (see Appendix A) in the 2D case. After eliminating 

v from the system (jSj) and going to non-dimensional variables (see Appendix B for the 
details) we obtain the following equation for q(t) 

? = -e-s+<7, (9) 

with e = — ^~~t i / t where p s is the density of the solid disk and q c har is the 

characteristic value of g. 

In the same way, we can consider the analogous three dimensional problem setting 
Q to be the half-space {(x,y,z), z > 0} and the particle to be a ball of radius R. The 
asymptotic expression is well known in this case (cf. Remark [TBI) and is given by 67rVu^-. 
Performing the same non-dimensionalizations as in the 2D case, we arrive at the equation 
for q(t) (the distance from the particle to the wall): 

q=-e-+g, (10) 
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Table 1: Typical values of the parameter e in © or in ([TO]) taking the viscosity and density 
of either water or glycerin and different value of the particle radius R and density p s . 

We remind that equations of the type (!9~])- (TTD"]) are at the basis of the gluey particle 
model of [H] [12]. The model consists in fact in considering the limit e — > 0, which is 
physically the limit of vanishing viscosity. We see, however, that e is not necessarily small. 

3.2 Estimates and an approximated solution for the model ODE 

Consider the ordinary differential equation 

q=~n(q)q + g (11) 

or, equivalently, the system of first-order equations 

v =- n (q)v + g 
q = v 
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where n(q) is a given differentiable decreasing positive function on (0, oo) with n(q) = 
N'(q) such that N(q) — > — oo as q — > 0. Note that this is valid for n(q) = s/q^ or 
n(q) = e/q, which give the asymptotic of the lubrication force in 2D and 3D respectively. 
We assume g G L ; 1 oc (IR + ) from now on. Under this hypothesis, Problem ffTTj) . completed 
with appropriate initial conditions, is well-posed, as proved for instance in [HJ Prop. 1.1] 
for the case n(q) = e/q. The proof is generalized without problem to any n(q) satisfying 
the above hypotheses, as mentioned in [HJ Section 5.1]. 

Proposition 1 Given (qo,Vo) G M 2 withqo > there exists a unique positive global solution 
q G W^(M, + ) to (fiTi) with initial conditions: 

q(0) = q q{0)=v o . (13) 

Let us take some threshold value q s . We suppose that q goes below this value after 
some time t\ when q(t) hits the threshold q s for the first time. Our aim is to find a suitable 
approximation of q after ti which enables to predict the time t 2 when q goes above the 
threshold q s without solving ( TTT|) . To this end, we introduce 

v(t) = qih) + f g(s)ds (14) 

for t > t\ and note that v{t 2 ) = qfo). Indeed, integrating ( fill from t% to t > t\ shows that 

q(t) = qih) + N(q s ) - N(q(t)) + f g(s)ds = N(q s ) - N(q(t)) + v(t). (15) 

Setting here t = t 2 and noting that g(t 2 ) = q s gives the desired result. Our approximation 
of the trajectory q(t), denoted by q(t), stems from the assumption (verified afterwards in 
Proposition [3]) that the velocity g(t 2 ) at the return point t 2 is small provided the threshold 
q s is small. Consequently, a good approximation for the time t 2 should be provided by the 
time i 2 defined as the first time larger than t\ when v(t) = 0. The construction of the 
approximated trajectory is hence the following: we first assume that q(t) is the same as 
q(t) until the latter hits q s for the first time at t = t\. Next, the trajectory q(t) is frozen 
until the time t = i 2 and resumes then again as a solution to ffTTl) starting from q s with 
zero velocity: 

(q(t), for < t < h 
q s , for ^ < t < i 2 (16) 
solution to ffTT]) with q(t 2 ) = q s , q(t 2 ) = 0, for t > t 2 

(see Figure El). 

In order to study the error committed by introducing the approximated trajectory q(t), 
we start by inferring the following estimate: 

Proposition 2 The time i 2 provides a lower bound for t 2 , i.e. i 2 <t 2 . 
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Figure 2: Construction of the approximated trajectories q(t) and q(t). 



Proof. There are three cases. If t 2 does not exist, there is nothing to prove. If t 2 = t\ this 
means t\ is a local minimum for q so that q(ti) = and t\ = t 2 = t 2 . Finally, if t 2 exists and 
satisfies t\ < t 2 , since q(ti) = q(t 2 ) = q s the function q(t) has a local minimum i min G (ti,t 2 ) 
where q{t min ) = 0. Evaluating ( 1T5|) at t — t min shows then v(t m i n ) = N(q(t m i n )) — N(q s ) < 0. 
On the other hand, v(t 2 ) = q(t 2 ) > since the function q(t) is increasing at t 2 by the 
definition of t 2 . The intermediate value theorem tells now that 3t G [t m m, ^2] such that 
v(t) = 0. Thus the first time t = i 2 when v(t) = is certainly in (t 1 ,t\ C (t 1 ,t 2 \. 
■ 

We have proved, in particular, that if the true trajectory q(t) ever returns to the 
values above the threshold q s , i.e. t 2 < 00, the indicator i 2 will show it, i.e. i 2 < 00. 
In Proposition El we compute an upper bound t 2 for t 2 . In particular, under suitable 
assumptions on g the time t 2 exists. 

To see moreover that i 2 may be a good approximation to t 2 we need, as mentioned 
before, to control q(t 2 ). The following proposition proves that q{t 2 ) is not too large with 
respect to q s . 

Proposition 3 Let q + (resp. g + ) be the positive part of q (resp. g): q + = max(g, 0) 
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(resp. g + = max(g , 0) ) . Then 



* (17) 



where \\ ■ \\ denotes the supremum on [ti,^]- particular, q(t 2 ) > ; so t/iat 5^2) — l^gT - 

Proof. Let q be non-negative on an interval [t m i n , t max ] C [ti, £2]- Without loss of generality, 
we can thus assume that £ m j n is a local minimum of q(t) and t max is either a local maximum 
or t max = t 2 . Then, we introduce 

/(t) := g(t) - J^L, VtGlti.ta]. 

If /(to) > at some to G (t m i n ,t max ) then /(t) rests positive on a small interval around to 
so that q < by invoking ( ITTj) and thus q(t) is decreasing on this interval. Since q(t) is 
increasing on [t min ,t max ] and n(q) is decreasing, we obtain that f(t) is decreasing locally 
around t so that f(t) > /(to) > for all t G [to — 5, t ] with some 5 > 0. Repeating this 
argument for times down to t m i n we see that /(t) should be positive on [t min , t ] but this is 
impossible since q(t m i n ) = and /(t m i n ) < 0. This proves that /(t) < on [t min , t max ]. ■ 
The results above allow us to fully characterize the error of approximating q(t) by q(t) 
up to time i 2 and to bound it by certain quantities depending on q alone for time after t~2, 
at least in the special case when g(t) is given by 



g(t) 



g-(t)<0, fort<t , 
g + > 0, for t > t , 



with some positive constants t , g+ and a negative function g~(t). Indeed, we prove in the 
following proposition that the error r = t 2 — i 2 is bounded by l/n(q s ) and we remind that 
n(q s ) — > 00 as q s — > 0. 

Proposition 4 If g(t) is given by (fT8j) . t/ien 

1 



< to - to < 



n{q B y 

sup \q(t)-q(t)\<q s , (19) 



te[0,t 2 ] 

0+ 



|g(t 2 )-^(t 2 )|=g(t 2 )< 



Moreover, both q(t) and q(t) are non-decreasing for t > i 2 and 

(q(t)-q(t-T), ift>t 2 , 

with t = to — to. 
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Proof. We first note that t 2 > t 2 > to- The first inequality here is already proved in 
Proposition [2j The second inequality follows from v(i 2 ) — 0, which can be rewritten as 



ft 2 

-Q(ti) = / g(s)ds. 
Ju 



Since q(t\) < 0, the last equality cannot hold if g is negative everywhere on [ti,^]- Since 
t 2 >t , we have g(t) = g + > for t >i 2 . Evaluating (fT5l) at t = t 2 gives 

q{t 2 ) = I* g(s)ds = g+(h-h). (21) 
Jt 2 

Recalling Proposition [3], we see that 

?(*2) / 1 



t 2 -t 2 = < 



9+ n(q a ) ' 



which is the first inequality in f|T9|) . The second inequality in f|T9|) is obvious since < 
q(i 2 ) < q s and q(i 2 ) = q s . The third inequality follows from Proposition [3] since q(i 2 ) > 0. 
Indeed, f lT5|) evaluated at t = i 2 gives ^(^2) = N(q s ) — N(q(i 2 )). 

We now turn to the study of q(t) and q(t) for time t > i 2 , in order to prove that 
q(t) and q(t) are non-decreasing on this interval. We first observe that q(t) cannot have 
local maximums at t > to. Indeed, at any local extremum t e > to, equation ( TTTf) would 
imply q(t e ) = g + > 0. Since q(t) is non-decreasing at i 2 > t , it should be non-decreasing 
also everywhere on [t 2 , 00). The same reasoning applies to q{t). The estimate (|20|) is now 
evident in the case i 2 < t < t 2 since q(t) > q s and q{t) < q s for such t. 

In order to prove (1201) in the case t > t 2 we will show that q(t) is squeezed between 
q T {t) and q(t) for t >t 2 where q T {t) = q(t — r). Indeed both q T (t) and q(t) satisfy the same 
equation f lTT]) with the same initial condition at t — t 2 (q(t 2 ) = q T {t 2 ) = q s ) but possibly 
different initial velocities q(t 2 ) > 0, q T (t 2 ) = 0. Writing ffTT|) for q T (t) and q(t), taking the 
difference thereof and integrating from t 2 to t yields 

q(t) - q T (t) = q(t 2 ) - N{q{t)) + N(q T (t)). 

This shows that if the trajectories of q(t) and q T (t) intersect at some time tl then qit!) > 
q T (t') so that q(t) > q T (t) at least for some time after t! . Since q(t) > q T (t) holds also for 
some time after t 2 it should hold for all t G [t 2 , 00). 

It remains to prove that q(t) < q(t) for all t G [t 2 , 00). To this end, we integrate ( ITT]) 
for q(t) from t 2 to t and that for q(t) from i 2 to t. This yields with the aid of ( f2~Tl) 



q{t) = q(t 2 ) + N(q s ) - N(q(t)) + [ g(s)ds = N(q s ) - N(q(t)) + I g(s)ds, (22) 

Jt 2 Jt 2 

q(t) = N(q s ) - N(q(t)) + [_ g(s)ds. (23) 

Jt 2 
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We see now that if q(t') = q(t') at some time t' > t 2 then also q{t!) = q(t'). It means by 
the uniqueness of solutions to (TTTj) that the trajectories of q(t) and q{t) either coincide or 
do not intersect on [i 2 , oo). Since q(t 2 ) < qii 2 ) this implies q(t) < q(t) on [i 2 , oo). 
■ 

It is possible to relax the hypotheses of the last Proposition on the particular form of 
the function g in several ways, if we still remain in the case when q goes below the threshold 
value q s a finite number of times (in particular, g should be allowed to change sign only a 
finite number of times). Although we do not have an analogue of these results for a general 
git), we can always provide an easily computable upper bound for t 2 alongside the lower 
bound t 2 . 

Proposition 5 Let t 2 be the first time after t 2 such that 

-t 2 

(f 2 - s)g(s)ds = q s . (24) 

t 2 

Then t 2 > t 2 . 

Proof. Consider q(t) defined for t > t 2 as the solution to 

q(t) = f g(s)ds, q(i 2 ) = 

We recall now the property (I2"2"j) of the original solution and note N(q(t)) < N(q s ) for 
t € {t 2 , t 2 ) so that 

q{t) > I g(s)ds > l(t). 

Since q{t 2 ) > qii 2 ), the trajectory of q(t) lies above that of q{t) for t G it 2) t 2 ) so that g(t) 
hits the threshold q s before g(t). In other words, t 2 < t 2 , where t 2 is the first time after i 2 
such that g(t) = q s . This definition of t 2 is equivalent to that in (124]) . ■ 



4 Discretization in time of the model ODE 

4.1 Three schemes for ODE (1111) 

We first describe a straightforward Euler discretization in time of (11 II) rewritten as the 
system (fT2l . Introducing the time step At and the discrete times tk = kAt, k = 1,2, .. . 
we thus consider the following 

Algorithm [1] is inspired by the real fluid-particle simulations in which one first finds 
the new velocity at each time step and then moves the particle with this velocity. An 
immediately evident drawback of the scheme in Algorithm [T] is that it does not necessarily 
provide a positive approximation q k . We remind that negative values of q are unphysical. 
Moreover, even if approximations q k remain positive but become tiny at some time steps, 
one would require very small stepsize to obtain an accurate solution. 
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Algorithm 1 (straightforward Euler) 



Step 0. Initialize (v°, q°). 

Step 1. For k = 1, 2, . . . update (v k , q k ) as follows 



= - n (q k -i)v k + g(t k ) 



At u ■ 



(25) 



Remark 6 One may argue that this can be cured by resorting to a fully implicit scheme 
that couples the evaluation of the velocity with the displacement of the particle: 

= -n(q k )v k + g(t k ) 



At 



(26) 

At ~ = V "- 



1 k -Q k 1 _ „* 



Indeed, this scheme provides a positive solution for q, for example, in the important case 
n(q) = e/q. To see this we eliminate v k from f )26|) . which gives an equation for q k only 

^- £ -^ = q -^-e + v k -i + Atg(t k ) 
At q k At yV ' 

The function of q k in the left-hand side is increasing form — oo to oo when q k goes from 
to oo. This means that there is the unique positive solution q k for any given v k ' 1 and 
q h ~ l . However, scheme ( 1261) would be too expensive in a real fluid/particle simulation where 
n(q k ) should be recomputed after any change in q k through a numerical approximation of 
the Stokes equations. 

Remark 7 Difficulties related to the discretization of these differential equations have al- 
ready been pointed out in [T^j (see Remark 3.1): it is shown that very small values of q can 
be reached (below the smallest real number which can be stored by standard numerical soft- 
wares and also below intermolecular distances). Different approaches have been proposed 
in previous related works to deal with similar problems: in Section 5.2], the author 
suggests the introduction of a cut-off function for the microscopic distance 7 = £ln(g), 
defined in terms of the roughness of the surface (see also \12\ Section 2.5]); in IT2\ Section 
2.4-3], for the case of fluid/particle simulation, a constraint for the distance q is defined in 
terms of the mesh size of the fluid domain. 

In order to overcome these difficulties, in this work we propose to introduce a small 
threshold value q s > 0, to replace the exact solution q(t) by the approximated one q(t), as 
summarized in ffTB"]) on the continuous level. Discretization in time of this idea introducing 
also an approximation of v(t) defined in (Til| is detailed in Algorithm [2J Note that the exact 
threshold is replaced in Step 2 of this algorithm by the last available value of q k before 
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Algorithm 2 (with an a priori fixed threshold) 



Step 0. Initialize (v°, q°) and choose q s > 0. 

Step 1. For A; = 1,2,... update {y k , q k ) by ( 1251) until q k < q s at some k = k\. We then 
redefine v kl = 0, q kl = <? initialize u by v kl = v kl ~ x + g(t kl )At and switch to 
Step 2. 

Step 2. For = ki + 1, fc x + 2, . . . update v k as follows 

tj* = v*- 1 + 5 (*fc)^* (27) 

keeping q k = q kl ~ 1 , v k = until w fc > at some k = k 2 . We then abandon the 
calculation of v and switch to Step 3. 

Step 3. For k = k 2 + 1, k 2 + 2, . . . update (v k , q k ) as in Step 1. If q k goes again under q s 
at some time step, then reintroduce v and keep switching between Steps 1 and 2 
back and forth. 



passing below q s , i.e. q kl 1 , so that the passage time t 1 is legitimately approximated by 
tki-i- 

As suggested by Proposition H] a good choice for q s would be such that l/n(q s ) = CAt 
with some constant C so that the criterion for switching from Step 1 to Step 2 can be 
rewritten as n(q k ) > l/(CAt). Indeed, extrapolating the result of Proposition H] to a 
general g{t), we see then that the error of approximation of the exact return time t 2 by 
i 2 ~ k 2 At is of order At. Moreover, the error \q(t) — q(t) \ is dominated by \q(t) — q(t — r)| 
with r of order At. Thus the error caused by the introduction of the threshold should 
be of of the same order as that of Euler scheme itself. However, an optimal choice of the 
constant C above would require some a posteriori error indicators to control the error of 
the Euler discretization. 

Generally speaking, a discretization of an ODE using the constant time step is often not 
optimal concerning the CPU time needed to achieve a desired accuracy. This is especially 
true for our ODE ffTTl) since the coefficient n(q(t)) can change enormously during a sim- 
ulation. One could try therefore to modify our algorithms by introducing a non constant 
stepsize chosen at each step using an error indicator. A general recipe for an optimal adap- 
tation of the stepsize, according to [lj is to choose the stepsizes so that the discretization 
error per time step remains approximately constant during the whole simulation. As is 
well known, the error per time step in a Euler scheme like that of Algorithm 1, is given 
to the leading order by max(|g(t fc )|, |ii(i fc )|)At^/2 where At k = t k — t^-i is the stepsize on 
step k which can be varying. Assuming that (q k ~ 1 , v k ~ l ) and (q k , v k ) are Euler approxima- 
tions at times tk-\ and t k respectively, the second derivatives of q and v can be computed 
approximately by deriving the equations in f|T2|) with respect to time and replacing the 
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Algorithm 3 (with adapted variable stepsize and automatically chosen threshold ) 



Step 0. Initialize (v°, q°) for to = 0, choose the first time step Ati > 0, the minimal 
tolerated time step At min > and the error tolerance per time step tol > 0. 



Step 1. For k = 1,2, . . . set t k = t k _i + At k , update (v , q ) by 

—n(q k ~ 1 )v k + g{t 



At 



fc-i h " ' ^ " (28) 



q k -q K ~ 1 _ „,fc 



and calculate the approximation for the error committed on this time step by 
e k = max(|g(£fc)|, |i)(£&)|)A^/2 where g(tfe), i5(ijfc) are computed using (12"91) - (13"01) . 
If e k < tol and > then we accept the just computed values {y k , q k ) and 



proceed to the next time step increasing the time step to At k+ \ = ^jiAt k . 



Otherwise, if e k > tol or q k < 0, we reject the approximation (v , q ) and try to 

recompute it by (|28|) with a smaller time step At k := min( ^J^-At k , At k /2). If the 

new approximation (v k , q k ) is still not sufficiently accurate, we try to diminish the 
time step again and again until an acceptable approximation (v k , q k ) is computed 
and proceed only then to the next time step taking At k+ i equal to the smallest 
value of At k used on step k. However, if in the process of reducing At k we come 
to a time step smaller than At m ; n at some k — ki, we abandon the calculation of 
(v kl , q kl ) and switch to Step 2 taking the initial values v kl ~ l = v kl ~ x , q s = q kl ~~ x 
and setting At to the current value of At k . 

Step 2. For k — k\, k\ + 1, . . . update v k by (|27j) keeping q k = q s , v k = until v k > at 
some k = k 2 . We then abandon the calculation of v and switch to Step 3. 

Step 3. For k = k 2 + 1, k 2 + 2, . . . update {y k , q k ) as in Step 1. If q k goes again under q s 
at some time step, then reintroduce v and keep switching between Steps 1 and 2 
back and forth. 



missing derivatives by finite differences. This gives 

?(«*)« -n(«V+j(« (29) 
«(*») » nW ~ n( a %[t k) - "(^-"(^V + »M -»<«>-,) (30) 

Thus, denoting E k = max(|g(t fc )|, \v{t k )\)/2 where the second derivatives are replaced by 
the formulas above, the strategy to choose the time step would be to require E k At\ = 
tol = const where tol is prescribed tolerance. As this recipe can lead to extremely small 
At k we combine it with the threshold approximation by switching to q(t) only when At k 
becomes smaller than some minimal admissible time step size. This idea is implemented 
in Algorithm [3j 



14 



4.2 Some numerical tests 



As an illustration, we consider the ODE ( fTTT) with n(q) = e/q 3 ^ 2 corresponding to the 
lubrication force in 2D. We report several numerical results setting 



g(t) 




for t < 2, 



, for t > 2, 



and the initial conditions g(0) = 1, q(0) = 0. In all the cases, we have at our disposal a 
very accurate reference solution, which we call "exact" in what follows. 

We first compare Algorithm 1 (without threshold) vs. Algorithm 2 with the threshold 
such that n(q s ) = 1/(20 At). The choice of the coefficient 20 in the last formula is somewhat 
arbitrary, but it works fine in our test cases. The first series of numerical experiments is 
performed taking e = 0.1. The results are reported in Figure [3j Note that, although the 
solution obtained with Algorithm 1 is positive at At = 0.01 and smaller, the introduction 
of the threshold in Algorithm 2 seems to enhance the quality of the solution at large time. 




Figure 3: Solution to the ODE ( TTTj) with n(q) = e/q 3 ^ 2 , e = 0.1. Left: the exact solution 
and solutions obtained by Algorithm 1 with At ranging from 0.1 to 0.0001. Right: solutions 
obtained by Algorithm 2 with At in the same range. The results obtained with At = 0.0001 
are visually indistinguishable from the exact solution. 



Let us now turn to another series of numerical experiments taking e = 10 -3 . The results 
are reported in Figure HI Algorithm 1 is now unable to produce a physically acceptable 
solution even with very small At = 0.0001. Algorithm 2, on the other hand, works fine. In 
order to observe better the quasi-contact region and the effect of introducing the threshold, 
we zoom on small distances by passing to a log-scale in Figure [5] (left). We report also some 
numerical experiments aiming at the determination of an optimal value of the threshold q s . 
As said already, our strategy is to choose it so that n(q s ) = l/(CAt) taking C = 20 in all 
the preceding experiments. Figure (left) illustrates the results obtained fixing At = 0.001 
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Figure 4: Solution to the ODE (QT]) with n(q) = e/q 3/2 , e = 10" 3 . Left: the exact solution 
and solutions obtained by Algorithm 1 with At ranging from 0.1 to 0.0001. Right: solutions 
obtained by Algorithm 2 with At in the same range. The results obtained with At < 0.01 
on the right are visually indistinguishable from the exact solution. 





Figure 5: Left: same results as on Figure HI i.e. varying At and choosing q s so that 
n(q B ) = l/(20At). Right: same test with At = 0.001 fixed and different values of q s 
calculated so that n(q s ) = l/(CAt) and C varying in the range [1,2000]. 

and taking three different values for q s by varying C. These results confirm two general 
observations: taking q s too small may deteriorate the accuracy of the solution (indeed, the 
original goal of the introduction of q s was to avoid the extremely small values of q). On 
the other hand, taking q s too large may unnecessarily perturb the solution in the regions 
where a straight discretization could give more precise results. Thus, there should be an 
optimal choice of q s which seems to be not too far from the formula mentioned above with 
C = 20. 



16 




Figure 6: Results obtained by Algorithm 3 with tol = 10~ 5 and different values of Ai m j n 
specified in the legend. Left: evolution of q(t k ) as compared to the exact solution. Right: 
evolution of Atk- 

We conclude by a series of experiments using Algorithm 3 with adapted stepsize. We 
take again the same governing equation with e = 10~ 3 and run Algorithm 3 setting the 
tolerance to tol = 10 -5 and taking a number of values for At m j n . The results are reported 
in Figure |6] with the evolution of q(t h ) on the left and that of Atk on the right. We give 
in particular the results with At min = so that Algorithm 3 is reduced to a standard time 
marching scheme with automatically adapted step sizes. It is interesting to note that step 
sizes At of the order 10~ 4 are actually sufficient to satisfy our error tolerance criterion 
almost everywhere apart from a very small range of time around 1 where the Algorithm 
chooses At of the order 10 -8 . Setting At m i n to 1CT 3 or 10 -4 avoid such small time steps 
and gives fairly good results. 
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A An estimate for the drag (lubrication) force 

This part is devoted to computing the first order expansion of the viscous drag exerted by a 
fluid on a disk approaching a plane wall. Let Q C M 2 be the half plane Q = {(x, y) with y > 
0}, the boundary of Q is dfl = {(x, 0), x G M} (see Figure H]). We assume that the fluid 
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fills the domain T q := ft \ B q where B q := B((0, R + q), R) is the domain occupied by the 
solid disk of radius R > and placed at the distance q > from the boundary is dQ. Our 
notations are lightly different here from that of Section [2l Indeed, domains T q and B q are 
now indexed by q rather than t since the eventual dependence of q on time does not interest 
us in this Appendix. Assume moreover that the disk has a prescribed velocity V = Ve2 
and zero angular velocity. The motion of the fluid is governed by Stokes equations (JT|) 
with g = and zero boundary conditions at the infinity for u. It is classical that (JT|) 
is well-posed in a suitable framework and has a unique solution (u,p) (the pressure p is 
defined up to a constant). In particular, the viscous drag force F exerted by a fluid flow 
(u,p) on B q and computed by ([2]) is well-defined as a function of V and q. Our aim is to 
prove the following result: 

Theorem 8 The viscous drag force satisfies F = — n(q)Y where 



with e(q) — > when q — >■ 0. 

Remark 9 The equivalent result in the three-dimensional setting is well-known but it is 
new in two dimensions. Whatever the novelty of this result, the main contribution in this 
section is that we provide a new way for computing the first order expansion of the viscous 
drag. This new method is more robust than the known computations |U GJ EES. HE/- In 
particular, we claim it extends to the three-dimensional setting with small changes and 
can be adapted to other boundary conditions. As an illustration, our method is a well- 
designed tool for estimating the distance between solid bodies in solutions to the fluid- 
structure interaction system with the full newtonian Navier Stokes equation on the fluid 
domain JB, El 13/ ■ 

The proof of Theorem [8] is divided into three steps. First, we recall the variational 
formulation for the solution (u,p) to (CEJ) and apply it to compute the associated drag 
coefficient n(q) in the expression for viscous drag F = —n(q)\. Secondly, we deduce from 
the variational formulation the lower and upper bounds for n(q). We conclude the proof 
by computing an asymptotic expansion of the bounds of this range. 

1. The variational formulation for the Stokes system ([TJ). As (u,p) in ([TJ) and 
the formula ([2]) for the drag depend linearly on V, we set V = e2 in what follows. It is 
classical to extend the velocity u to the whole domain Q by setting it to e2 inside B q and 
to look for u in the following space 



The solution of ([T]) is associated with the minimization of Dirichlet integral D(w) = 
J n |V(w)(a;, y)\ 2 dxdy over the subset Y q of Dq 2 (Q): 





w G Ll oc (Q) with Vw G L 2 (tt), div w = in and w = on dtt 
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Let us denote this minimum by Ed(q), i-e. 



E D (q) := min {D(w) ; w G Y q } , (31) 

As the Dirichlet integral D is a strictly convex functional on Y q it has a unique minimizer 
u g for any q > 0. It is easy to see that u q gives the solution of (JT]). Indeed, we have for 
any w G C^°(J-q) such that div w = 0, 



/ Vu(i, y) : Vw(x, y)dxdy = 0, 
in 



which is the variational formulation of (JTJ. The pressure p q G Lf oc (Q) can be then recovered 
as the Lagrange multiplier corresponding to the constraint div u = 0. We notice that p q 
is unique up to a constant and the ellipticity of the Stokes problem implies that (u q ,p q ) 
is smooth on J- q so that it furnishes a classical solution to ([T|). We refer to [I] for more 
details and also for a proof of the converse implication. 

Given w G C^°(Q) such that div w = and w = e 2 on B q , a straightforward integration 
by parts yields (recall that n is the unit normal looking into J r q ): 

v I Vu q (x,y) : Vw(i, y)dxdy = 2v \ D(u q )(x,y) : Vw(x, y)dxdy 

= - [2uD(u q ) - p q I 2 }nda ■ e 2 . 

JdB q 

Taking a suitable family of approximation of u q (see jl] for details) we obtain in the limit: 

vE D (q) = - / [2vD(u q ) - p q I 2 ]nda ■ e 2 . 

JdB q 

The drag F is parallel to the velocity V for symmetry reasons, i.e. F = (F • e 2 )e 2 . We 
have thus the following result. 

Lemma 10 For any given q > and V parallel to e 2; the drag force is given by F = 
— n(q)Y with n(q) = vEo^q). 



2. Upper and lower bounds for Erj(q). We apply the above variational formulation 
to bound En{q). To this end, given q > 0, we first prove that the Dirichlet integral of any 
u G Y q is greater than some mn^q) depending only on q. This gives us a lower bound for 
Er>(q). Then, we construct a suitable u q G Y q and compute its Dirichlet integral Mc(q). 
This gives an upper bound for Ejj(q). We finally compare and Mjj for small values of 
q to prove Theorem [81 

As the scale invariance of our problem implies that n(-) is actually a function of q/R, 
it is sufficient to consider the case R = 1, which we admit in what follows. For any 
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x G (—1, 1), y < 1 + q such that (x, y) G dB q there holds y = 1 + q — Vl — x 2 := 5 q (x) and 
we denote 

tt q := {(x,y) G M 2 such that x G (-1/2,1/2) and y G (0,<L(a;))}. 



Given a smooth u G Y q , we introduce 

^(x,y) = 



ui(x, z)dz. 



Then, straightforward computations imply that ip G C 00 (J r g ) and satisfy u = V -1 ^ 
(—dyip,d x ilj), so that we have the boundary conditions: 

f V^y) = 0, d y 4>(x,y) = 0, onffl, 
\ d x ip(x,y) = 1, d y ip(x,y) = 0, on dB q . 

On the other hand, we compute the Dirichlet integral D(u) with respect to ip. This yields: 



(32) 



u > 



d yy il)(x,y)\ 2 dxdy. 



We denote I (if}) the integral on the right-hand side of the above inequality. Thus, D(u) is 
greater than the minimum of I (if}) over smooth ip satisfying the boundary conditions ( 132|) . 
This minimum is computed in the following lemma: 

Lemma 11 Given q > andip G C°°(Q q ), satisfying boundary conditions (1321) there holds 
(VO > H^g) w here : 

2 



i/) q (x,y) = x 



y 



_5 q (x) 



y 



AO). 



, \/(x,y)ett q . 



Proof. Given any ip(x,y) G C°°(Q q ), satisfying ( l32|) . boundary conditions on dB q imply 
if>(x,5 q (x)) = C+x,\fx G (—1/2, 1/2), with some C G R. Thus i\) satisfies Vx G (-1/2,1/2) 

^(x,0) = 0, V(^^)) = C + ac, ^(x,0)=0, ^(x, 5 g (x)) = 0. 

Then, for arbitrary a; G (—1/2, 1/2), excepting eventually a; = — C, we introduce Xx(t) = 
ip(x,t5 q (x))/(x + C), t G (0, 1), which satisfies 

X,(0) = 0, Xx(l) = l, Xi(0) = 0, x / *(l)=0. (33) 

By a standard optimization argument we see that the minimum of \x'x(t)\ 2 dt over a ^ 
smooth Xx(t) satisfying ( 1331) is attained on the function r](t) such that r]^'(t) = on (0, 1), 
which is given by rj(t) = t 2 (3 — 2t). This yields 

^ 2 ^(x + C) 2 



\d yy ip(x,y)\ dxdy 



> 



-1/2 JO 
1/2 /•! 



1/2 JO 



{x + Cf 

IWI 3 



WL{t)\ 2 dtdx 



\rj"(t)\ 2 dtdx. 



(34) 
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Minimizing the last integral with respect to C, we obtain C = 0. The last inequality 
becomes equality if we take if)(x,y) = ip q (x,y) = xrj(y/5 q (x)). ■ 
Finally, a lower bound is given by 



m D (q) = \dyyip q (x,y)\dxdy. 

On the other hand, it is always possible to extend ip q on the whole Q q to some ip q G 
C 00 (J-" g ) nC°°(B q ) nC(fl) and such that u q := Vip q G Y q . For instance one might interpolate 
(in the x-variable) ip q with a suitable truncation of ipo(x,y) = x. As, outside Q q , the solid 
disk remains at a positive distance of dQ, the truncation and interpolation function can 
be made independent of q. In particular there holds: 

Vu q (x,y)\ 2 dxdy = m D (q) + r(q) + / [\d xx ip q (x, y)\ 2 + \d xy tp q (x, y)\ 2 ] dxdy := M D (q), 
n Jn q 

with r a bounded function of q. 

3. Asymptotic expansion of m£>(q) and Mo{q)- 
Lemma 12 There holds, for small values of q : 

m D (q) = 3v^ri±iM, (35) 

with e(q) —> when q — > 0. 

Proof. As already seen in the proof of the preceding Lemma, m£>(q) is given by the integral 
( |34|) with C = and = t 2 (3 — 2t). Integrating with respect to t and performing a 
change of variables x — >■ y/gx yields 

f 1 / 2 , 12 Z" 1 /^) x 2 dx 

m D (q) = 12 / Mq dx 



'-1/2 l<^0)| 3 " <? 3/2 J-i/M {5 q {^/qx)/qy 
Expanding to the first order, we have: 
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lim — — , , , „ = lim 



i^{5 q {y/qx)/qy 9 -o [(1 + q _ _ gx 2 )/g ]3 (l + x 2 /2) 3 ' 

We remark moreover that standard geometric arguments imply that S q (x) > q + x 2 /2 or 
dq(s/q~ x )/q > 1 + x 2 /2 so that we can apply the Lebesgue theorem to obtain: 

lim f 1/M x2dx - r x2dx - 

■ 

One can also check that c^V^ is the second derivative of which diverges the fastest 
when q goes to 0. Thus, M D (q) — mn(q) = o(m^(g)) as q — > so that (j3"o"|) holds also for 
Mo(q), which proves Theorem [HJ 
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Remark 13 The same method can be easily adapted to the three dimensional case. More 
precisely, the result of Lemma [721 is still valid : the viscous drag exerted by a fluid on a 
sphere can be computed as the result of a minimization problem of the same type as the one 
given by equation (3l\) . Next, since the problem is axisymmetric, one can use cylindrical 
coordinates and perform the same asymptotical analysis as in the 2D case, in order to 
obtain that 

d2 

n(q) = 6ttVu— [1 + e(q)] , 

q 

with e(q) — > when q — > 0. 

Numerical evaluation of the drag force. 

In order to compare the theoretical estimate given by Proposition [8] to a numerical eval- 
uation of the force F, we carry out the following simulations: we begin with solving the 
Stokes problem in a fluid with viscosity v = 1 and then compute the corresponding force 
exerted by the fluid on a particle of radius R = 0.1 situated above a plane at a distance 
q. The computations are performed by using the Finite-Element solver FreeFem++ [7\ 
and the results are plotted in Figure (dashed line with squares). They agree with the 

3 

asymptotic expansion F ~ "i\pli\Vv C^j 2 (solid line). 

Lubrication force - 2D 
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Figure 7: Approximation of the lubrication force 



Lubrication force : asymptotic expression 

Lubrication force : numerical result —a- — 
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B Non-dimensional form of the governing equations 
(ISD involving the lubrication force 



Let us first consider the 2D case. Eliminating the velocity v from (JHJ) and specifying 
m = p s 7rR 2 with p s the density of the solid, we arrive at 

p s nR 2 q = -3V2ttu (~J Q + {Ps ~ P/K-^V 

In order to understand the typical values of parameters in this ODE, we should pass to 
non-dimensional variables, which can be introduced as follows: 

q = Rq, t = Tt', g = g cha rg' 

where the radius R is used as the length scale, the time scale is denoted by T and g c har 
is a typical value of the external force density, so that g c har ~ 10m/ s 2 . We choose then 
the time scale T so that the non-dimensional external force becomes of order 1, i.e. 
IT PS p Pf 9 char = 1 an d obtain the non-dimensional ODE (dropping the primes) 

The asymptotic in the 3D case is is given by n(q) = 67rv— (cf. Remark [T3|) so that 
starting from the dimensional ODE 

Ps ^7TR 3 q = -Q-kvR 2 ^- + {p s - p f )\-KR z g (3D). 



and using the similar non-dimensionalizations as before we obtain ffTOj) . 
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